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RAPID  GENERATION  OF  SYNTHETIC  SEISMOGRAMS  IN  LAYERED 
MEDIA  BY  VECTORIZATION  OF  THE  ALGORITHM 

By  Robert  A.  Phinney,  Robert  I.  Odom*,  and  Gerard  J.  Fryer 

The  theory  and  practice  of  computing  synthetic  seismograms  in  layered  media 
are  now  quite  well  understood*  and  are  nicely  set  forth  in  the  monograph  by  Kennett 
(1983)  and  the  textbook  by  Aki  and  Richards  (1980).  In  applications,  these  methods 
can  still  be  unsatisfactorily  slow. 

It  may  be  necessary  to  work  with  models  having  a  large  number  of  microlayers. 
Modeling  of  reflection  seismology  data  in  typical  sedimentary  rocks  may  commonly 
require  more  than  100  layers.  Many  situations  having  a  smooth  velocity  gradient 
and  need  many  microlayers  to  model  the  gradient,  when,  for  one  reason  or  another 
asymptotic  methods  are  not  suitable.  For  example,  if  we  are  studying  wide-angle  or 
long-range  prot>agation  in  the  top  500  m  of  a  sedimentary  section,  we  need  to  model 
the  pronounced  near-subbottom  material  gradients  due  to  compaction  as  well  as 
particular  layers  which  are  high  or  low  velocity,  density,  or  Q  anomalies  with  respect 
to  the  reference  gradient. 

In  seeking  to  determine  the  variation  of  material  properties  in  a  sedimentary  pile, 
we  may  start  with  some  generalized  reflection  data  set,  containing  reflected  and 
refracted  signals  produced  by  an  impulsive  point  source  and  recorded  at  a  number 
of  sensors  in  the  water  layer  at  various  offsets  from  the  source.  Any  true  inversion 
procedure,  regardless  of  details,  must  compare  the  observed  data  with  synthetic 
data  for  a  sequence  of  models  in  order  to  obtain  and  characterize  models  which  fit 
the  data  within  some  prescribed  limits.  We  thus  seek  a  forward  modeling  procedure 
for  arbitrary  layered  models  which  is  fast  enough  that  many  repetitions  are  possible 
in  a  reasonable  time. 

The  procedure  that  we  describe  in  this  paper  is  based  on  a  reorganization  of  the 
inside  loops  of  a  conventional  reflectivity  algorithm  to  permit  vectorization.  With 
this  optimization,  it  is  unncessary  to  be  concerned  in  detail  with  optimization  ot 
the  code  that  generates  individual  matrix  elements  from  their  algebraic  expressions, 
a  technique  adopted  by  Kind  (1976).  The  vectorized  procedure  is  by  itself  faster 
than  previously  described  procedures,  simply  because  practically  all  redundancy  has 
been  eliminated.  When  implemented  on  an  array  processor,  however,  it  achieves 
additional  speed  enhancement  of  about  20,  by  making  effective  use  of  the  vector 
hardware.  We  have  implemented  this  procedure  on  a  Data  General  Nova  3,  and  on 
a  Perkin-Elmer  3230,  with  a  CSPI  MAP  300  array  processor,  in  both  Fortran  and 
C.  In  particular,  we  have  used  a  slowness  method,  with  the  Haskell-Dunkin  6x6 
minor  matrices  for  construction  of  the  reflectivity  (Haskell,  1953;  Dunkin,  1965). 
An  optimal  code  in  terms  of  stability  and  computational  overhead  would  probably 
use  the  Kennett  2x2  reflection-transmission  matrices.  We  will  take  as  a  basis  for 
most  of  this  discussion  the  monograph  by  Kennett  (1983)  and  not  repeat  much 
standard  material  here. 

A  vectorized  version  of  the  WKBJ  method  has  been  developed  and  tested  by 
Chapman  and  Orcutt  (1985),  and  applied  to  the  inversion  of  marine  refraction  data 
by  Shaw  and  Orcutt  (1985).  For  certain  problems  requiring  summation  6f  many 
rays  to  synthesize  the  signal,  or  when  investigating  phases  omitted  by  the'  WKBJ 
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theory,  the  more  accurate,  but  also  the  more  time-consuming  reflectivity  method 
would  be  required.  This  has  been  a  motivation  for  developing  our  algorithm. 

The  Procedure 

Consider  an  impulsive  point  pressure  source  in  a  layered  stack  of  sediments,  with 
a  surface  water  layer,  and  the  observed  acoustic  pressure  field  near  the  surface.  In 
the  slowness/reflectivity  method,  the  acoustic  pressure  field  P(r,  t)  may  be  repre¬ 
sented  by  the  double  integral 


P(r,  t) 


due^'Slu)  I  dp  pJ0(wpr)${wp,  p) 


where  p(=  k/u>)  is  slowness,  4>  is  the  plane  wave  response  function  of  the  layered 
medium,  and  S(u)  is  the  Fourier  transform  of  the  source-time  function  5(f). 
Carrying  out  the  integration  over  to  in  (1)  gives  the  exact  relation  for  the  solution 
as 

P(r,  t)  -  \  S(t )  •  D2  fo  *(p,  r)  .  J0 r)  dp,  (2) 

where  J0(t)  is  the  inverse  Fourier  transform  of  </0(w),  and  i  is  the  inverse  Fourier 
transform 


’■r,-£L 


e~m‘S(u){’(p,  w)Rpp(up,  u)a(p,  w)  du. 


Although  we  consider  cylindrical  waves,  may  be  regarded  as  the  plane  wave 
transient  response  of  the  medium  with  the  given  source  and  receiver  functions.  To 
generate  the  plane  wave  transient  response,  we  need  the  following  components  from 


S{w)  =*  Fourier  transform  of  the  source  ,:v"  wavelet 
f(p,  w)  =»  source  directivity  function  (=  1/.  ^ in ,  f'  r  a  point  pressure  source) 
Rpp(wp,  <•>)  *  reflectivity  or  total  medium  plane  response  function 
a(p,  w)  «*  receiver  directivity  function  (=  1,  for  a  pressure  receiver). 

In  this  paper,  we  present  a  scheme  by  which  the  computation  of  the  plane  wave 
transient  response  i  ( p,  r )  may  be  optimized  for  use  with  an  array  processor.  The 
evaluation  of  (2)  is  a  separate  question,  for  which  no  particularly  rapid  algorithm 
has  been  proposed,  in  the  general  case.  In  the  very  common  instance,  however, 
where  the  source  receiver  offset  is  not  close  to  zero,  (2)  can  be  rapidly  evaluated  by 
i  he  slant  stack 


P(r,  t) 


±£  H*  r 

2  x  dt  Jo 


i(p,  t  +  px)  op 


(Chapman,  1978;  Phinnev  of  al.,  1981). 

In  order  to  clarify  our  new  algorithm,  we  summarize  the  familiar  spectral  and 
slowness  methods  {(A)  and  ( B ) ). 
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(A):  for  (w  ==  0  to  w  =  wm„) 


initialize 

for  (fc  =  0  to  k  =  ) 

( a)  fixed ) 

for  (layer  =  0  to  layer  -  N) 

(k,  a;  /ixed) 

evaluate  layer  and  interface  matrices 
matrix  multiplies 

integrate  over  k  (Bessel  function  evaluation) 
integrate  over  u>. 

Similarly,  for  the  sbwness  method  we  have 

<B):  for  (p  =  0  to  p  =  p^) 


initialize 

for  (w  «  0  to  m  *  m— ..) 

<p  fixed) 

for  (layer  *  0  te  layer  *  S) 

<  w,  p  fixed) 

evaluate  layer  and  interface  matrices 
matrix  multiplies 

integrate  over  w  ( Fast  Fourier  Transform  — ► <j>  ( p,  r ) ) 
integrate  over  p  ( slant  stack ) 

( B )  is  more  efficient  than  ( A  )  as  it  avoids  the  Bessel  function  evaluations  and  the 
oscillatory  integral  over  k.  Both  algorithms,  however,  have  essentially  the  same 
structure  for  the  innermost  loops.  Starting  from  the  slowness  algorithm  (B),  we 
recognize  the  inner  loops  to  achieve  a  vectorized  algorithm,  which  is  an  optimum 
way  of  generating  the  plane  wave  transient  response  4>(p,  t). 

The  principle  computational  burden  in  evaluating  <Mp,  r)  lies  in  the  evaluation 
of  the  response  of  the  medium  to  plane  waves  paramaterized  by  k  =  wp  and  <*>.  We 
are  required  to  evaluate  the  reflectivity  of  the  stack  of  layers,  as  seen  at  a  reference 
level  in  layer  0.  In  the  Haskell-Dunkin  minor  matrix  method,  which  we  have 
implemented,  and  which  we  discuss  for  concreteness,  a  six-element  vector 


Vy  =  col[l,  0,  0,  0,  0,  0) 


is  established  in  the  lower  half-space,  and  analytically  continued  upward  to  layer  0 
by  a  sequence  of  matrix  multiplications 


Wo  —  QoQl  •  ■  ■  Q\- lVy.  15) 

(The  snecific  form  taken  by  these  matrices  depends  on  the  exact  starting  definitions 
used  in  setting  up  the  problem.  See  the  Appendix  for  the  formulation  us$d  here.) 
\Ve  then  obtmn  directly  vhc  reflection  iuiicuouj  uv  inspection  of 


Vq  —  col(A,  Rps<\,  — Rss A,  Rpp A,  fl.s/>A,  det  /?A]. 


(6) 
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Since  the  reflection  functions  arise  from  u0  as  ratios  of  elements,  the  procedure  (5) 
may  be  modified  by  multiplying  vn  by  any  complex  scaling  constant  at  any  stage  n, 
without  affecting  the  result.  Each  wave  propagator  matrix  <j»„  is  the  product  of  an 
interface  matrix  Fn  and  a  layer  crossing  matrix  E„ 

Qn  =  EnFn. 

We  will  obtain  a  fast  algorithm  by  exploiting  two  properties  of  these  matrices 
F„  can  be  written :  F„(p) 

E„  can  be  written  (7a) 

diag(e'‘*^(,,"+',•  1,  1,  e (7b) 

where  d„  is  layer  thickness  and  qn,  q„ '  are  P  and  S  vertical  slownesses. 

We  can  exploit  the  properties  (7a)  and  (7b)  to  convert  (B)  to  a  vectorizable 
algorithm,  as  follows. 

1.  Loops  are  nested  (outer  to  inner):  loop  inp;  loop  on  layer  (n  J;  loop  on  frequency 

[m]. 

2.  The  innermost  ( frequency)  loop  is  redefined  to  consist  of  vector  processes,  in 
which  the  dimension  of  the  relevant  complex  vectors  is  the  number  of  frequen¬ 
cies,  Af,  being  evaluated.  On  an  array  processor,  this  loop  can  be  implemented 
as  true  vector  operations,  while  on  an  ordinary  processor,  the  vectorized  code 
provides  substantial  computing  advantages  over  the  conventional  procedure. 

3.  Suppose  that  the  code  is  to  be  run  for  M  equally  spaced  frequencies.  The  six- 
element  complex  state  vector  v  is  assigned  storage  for  a  6  x  Af  complex  vector 
v 

1st  frequency 
2nd  frequency 


Mth  frequency  (8) 

This  vector  is,  for  each  new  value  of  p,  initialized  with  =  1.0,  for  m  *  1, 
2,  •  •  ■ ,  Af,  and  all  ocher  elements  set  equal  to  zero.  The  calculations  then 
proceed  by  carrying  out  each  successive  matrix  multiplication  in  (5)  for  all 
frequencies  by  operation  on  the  extended  6  x  Af  state  vector  v. 

4.  Consider  the  matrix  multiplication 

vn  =  E„Fnv„+l  (9) 

in  which  the  state  vector  is  analytically  continued  from  the  top  of  the  (n  4-  1) 
layer  to  the  top  of  the  (n)  layer.  When  we  extend  u„  to  t>„,  by  including  all 
freouencies,  we  need  suitable  extensions  of  Fn  and  En.  Now,  the  complex  6  x 
6  matnx  F  is  (7a)  independent  of  frequency.  We  may  thus  provide  stojage  for 
the  values  of  F  for  the  N  interfaces,  which  may  be  precomputed  once  for  each 
value  of  p  (the  outer  loop).  The  matrix  multiplication 


vtT)  = 


ii  <i>  n  ui  in 

,  Vi  .  •  •  • ,  i 

ii  (2)  n  ,2>  n  ,2> 

ui  i  uj  ,  •  •  • ,  , 


„  (AO  „  (Ml 
fl  I  ^  I 


l>9 


(AO 
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uses  the  same  36  values  of  Fn  in  operating  on  0„  for  all  frequencies.  In  an 
array  processor  environment,  this  operation  can  be  set  up  as  a  complex 
convolution  of  the  rows  of  Fn  with  un,  with  a  shift  of  6. 

5.  Matrix  multiplication  by  E„  requires  that  E„  be  given  extended  storage  for  M 
frequencies.  We  treat  E„  as  a  complex  six  vector,  since  it  is  diagonal,  and 
obtain  a  6  x  Af  representation  similar  to  (8),  we  we  call  En.  The  matrix 
multiplication 


On  -  tnU  (10) 

becomes  a  dot  product  of  two  complex  6xAf  dimensional  vectors. 

6.  Computation  of  the  values  of  £n  requires  some  optimization,  for  the  algorithm 
would  seem  to  demand  6  x  M  x  N  complex  values  to  be  computed  and  stored. 
In  particular,  we  seek  to  reduce  what  appears  to  be  2  x  M  x  N  complex 
exponential  evaluations.  We  do  this  by  generating  the  values  for  En  recursively, 
using  the  first  six  elements  (zero  frequency)  as  starting  values.  Assign  storage 
for  N  complex  “incrementor”  vectors  Gn.  If 

(11) 

and  similarly  for  sn~,  then  we  define 

Gn  -  (l/s„+,  1,  l/a„-,  s„-,  1,  sn+].  (12) 

The  collection  of  Gn  depends  on  p,  but  not  on  w,  and  may  be  precomputed 
once  for  each  value  of  p,  in  the  outer  loop.  Consequently,  £n  is  generated  by 
assigning  starting  values  (for  zero  frequency):  E\,z  «  *  1.0,  and  by  using  G„ 
to  generate  the  higher  frequency  values  of  En  by  complex  multiplication 


"  Em' 

Em+l 

V 

Em~ 6 

E, n-5 

- 

• 

• 

G , 

Em- 1 

for  m  -  7,  13,  •  •  •,  6  X  (M  -  1)  +  1.  (13) 


The  basic  computational  unit  thus  consists  of  the  following  vector  process  (which 
we  label  ( P) )  which  has  the  effect  of  continuing  the  system  vector  0  from  the  top 
of  layer  n  +  1  to  the  top  of  the  layer  n,  for  all  frequencies 

•  Complex  convolve  the  rows  of  Fn  with  £>„+(,  with  a  skip  of  6  — »  u. 

•  Complex  dot  product  of  En  with  u  — »  o„. 

•  Generate  new  vector  £n-i  for  the  next  layer. 

These  may  be  implemented  as  array  processor  calls  or,  lacking  an  array  processor, 
as  an  inner  loop  over  frequency.  On,  0N-i,  •  •  -  ,  may  all  occupy  in  succession  the 
same  storage  array. 

The  algorithm  is  summarized  as  ( C ) 


<C):  for  (p  =  po  to  p  -  p^) 
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( compute  the  reflectivity  for  this  p,  all  w ) 

initialize  £ 

initialize  Gn,  for  all  n 

for  ( layer  »  N  -  1  to  layer  =  0) 

(continue  v  from  one  layer  to  next) 
execute  vector  process  ( P ) 

(handles  all  frequencies ) 

Generate  plane  wave  seismogram  by  Fast  Fourier  Transform. 

One  pass  through  the  stack  of  layers  produces  the  complex  reflectivity  for  all 
frequencies;  the  plane  wave  seismogram  is  then  produced  by  Fourier  transformation. 

Incorporation  of  source  and  receiver  directivity  and  spectrum  into  the  algorithm 
is  easily  done,  as  is  a  frequency  domain  operator  embodying  H*(d/dt)  [equation 
(4) ].  The  algorithm  may  be  easily  redesigned  to  use  one  of  the  more  popular  matrix 
schemes,  such  as  the  reflection-transmission  matrix  scheme  of  Kennett  (1981, 
1983),  which  is  wholly  free  of  potential  overflow  or  underflow  problems.  Incorpo¬ 
ration  of  the  effects  of  a  free  surface,  and  of  source  and  receiver  at  different  levels, 
requires  the  construction  of  a  response  matrix  from  the  reflectivity  matrix  using 
methods  described  by  Kennett  (1983),  but  lie  outside  the  inner  loop  and  do  not 
appreciably  affect  the  computational  load.  In  the  Harkell-Dunkin  method,  the 
computation  of  the  elements  of  £  can  lead  to  overflow  and  underflow  at  higher 
frequencies,  particularly  when  the  compressional  wave  function  in  one  or  more 
layers  is  evanescent.  To  avoid  these  problems,  we  normalize  £  by  its  largest  value; 
a  process  which  does  not  affect  the  reflectivities,  as  we  have  remarked. 

Speed  and  storage.  The  performance  characteristics  of  our  implementation  are 
summarized  in  Table  1.  Even  though  the  algorithm  is  based  on  formation  of  long 

TABLE  1 

Summary  or  the  Performance  Characteristics  or  our 
Implementation  or  the  Rznscnvn-y  Algorithm 

_ SpjWMid  Sto^ei  for  Alfonhm  (C) _ 

Processor  Perkin-Elmer  3230  with  lk  cache,  3 mb  memory, 

and  floating  point  hardware 

Array  processor  CSPI  MAP  300 

Operating  system:  Unix  version  VU 

Language;  C 

Teat  speed:  50,000  indexed  floating  point  multipiy-adds  per 

second 

Algorithm:  Tested  for  256  frequencies,  one  value  of  p 

2.2  tec  per  layer  3230,  with  subscripted  array 
references 

0.6  sac  per  layer  3230,  with  pointers  for  array 

access 

30  msec  per  layer  array  processor 

Storage:  Given  N  layers,  M  frequencies,  and  using  N  -  25, 

M  -  256. 


Array 

Hz 

Worth 

522 

V 

Complex 

6  x  M 

48  x  M 

12,288 

f. 

Complex 

36  x  N 

144  X  N 

3,600 

6 

Complex 

6  x  M 

48  x  M 

12,288 

G. 

Complex 

6  x  N 

48  x  N 

600 

Total 


28,776 
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vectors,  from  Table  1  we  see  that  the  storage  actually  required  is  still  well  within 
the  capability  of  commonly  used  computers.  Use  of  pointer  access  methods,  which 
are  available  in  the  C  language  (or  in  assembler  patches  under  Fortran),  provided  a 
speed  improvement  factor  of  3.6.  The  extra  factor  of  20  obtained  by  use  of  the  array 
processor  is  clearly  a  major  factor  in  making  this  algorithm  one  that  might  be  used 
in  repetitive  situations.  While  any  code  may  be  marginally  improved  by  close 
attention  to  eliminating  unnecessary  operations,  such  as  multiplication  by  1  or  0, 
this  procedure  appears  to  be  nearly  optimum.  Function  evaluations  (exponentials 
and  square  roots)  have  only  a  trivial  effect  on  the  total  time,  and  matrix  element 
evaluations  have  been  removed  from  the  innermost  loop,  as  far  as  possible.  We 
judge  that  any  further  order  of  magnitude  improvement  in  performance  will  be 
through  the  use  of  suitable  hardware,  such  as  a  state-of-the-art  (250  mflops)  array 
processor. 

Use  in  an  inversion  process.  We  seek  particularly  high  speeds  to  facilitate  the 
inversion  of  reflection-refraction  data.  In  an  inversion  procedure,  one  is  asked  to 
take  a  single  data  set  and  to  generate  a  potentially  large  number  of  synthetic  data 
sets  for  iteration  to  a  solution  and  for  assessment  of  the  class  of  models  which  fit 
the  data  in  some  sense.  An  optimum  data  set  for  this  purpose  would  consist  of 
sufficiently  many  receivers  in  a  radial  array,  so  that  spatial  aliasing  and  windowing 
effects  are  unimportant.  It  would  then  be  possible  to  generate  an  accurate  represen¬ 
tation  of  the  pla^c  wave  decomposition  of  the  data.  The  inversion  process  would 
then  generate  synthetic  data  for  a  very  restricted  subset  of  slownesses,  which 
contains  most  of  the  independent  information  required  to  fit  the  data.  Suppose,  for 
example,  that  five  plane  wave  seismograms  are  being  generated  per  pass  at  a  model, 
and  2  sec  are  required  per  slowness;  then  360  models  can  be  evaluated  per  hour.  If 
it  is  desired  to  generate  the  spatial  representation  4>(r,  t ),  then  over  200  plane  wave 
seismograms  are  needed  as  input  into  a  reasonably  accurate  slant  stack,  and  only 
nine  models  per  hour  can  be  generated. 

Discussion 

Our  purpose  in  needing  a  fast  algorithm  is  to  be  able  to  explore  a  large  variety  of 
models,  both  globally  and  incrementally.  At  this  point,  we  discuss  these  models  in 
the  sense  that  we  are  simulating  a  laboratory  experiment  and  trying  to  develop  a 
more  thorough  understanding  of  the  wave  propagation  phenomena  which  charac¬ 
terize  the  sea-bottom  sediment  column.  This  understanding  is  a  prerequisite  to 
applying  the  modeling  procedure  in  an  inversion  process,  which  will  require  judicious 
choice  of  procedures  for  subselecting  data  and  model  parameters. 

At  this  point  we  have  no  idea  whether  the  speed  achieved  is  fast  enough.  With 
cne  array  processor,  computation  speeds  which  deliver  tens  of  models  per  hour  are 
certainly  insufficient  for  even  a  restricted  Monte  Carlo  procedure;  perhaps  more 
orderly  inversion  algorithms  can  be  established  which  can  make  good  use  of  this 
capability.  With  modem  supercomputers  and/or  modem  superarray  processors, 
however,  hundreds  of  models  per  hour  are  now  feasible.  Whether  the  application 
justifies  that  kind  of  expense  remains  to  be  seen. 
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Appendix 

The  interface  matrix  Fn  is  the  matrix  product 

Fn  -  Tn~lT^t 

where  T„  is  the  6x6  delta  matrix  for  the  nth  layer.  The  elements  of  the  Dunkin 
matrix  are  formed  from  the  2x2  subdeterminants  of  the  4X4  matrix  which 
generates  the  stress-displacement  vector  irom  the  wave  potential  vector.  The 
interface  matrix  Fn  is  explicitly  frequency- independent,  and  because  we  are  inter¬ 
ested  only  in  ratios,  we  can  ignore  any  explicit  normalization,  and  common  factors 
among  the  elements  may  be  removed. 

The  following  is  a  list  of  the  16  independent  elements  pf  the  Dunkin  matrices. 
To  reduce  the  number  of  multiplications  required  to  form  the  elements,  they  have 
all  been  divided  by  p. 

tu  =  ~(p2  +  qq')/n  =  tie 

t12  =  —2pq/p 

tn  =  ~(p?  ~  qq')/p  =  -tu 

f15  =  —2pq'/p 

<2i  *  iq' IS2  —  —  to  —  —tu  —  —tie 
«3i  =  ~ip(  T  +  2qq  ' )  =  fa*  =  tu  « 
f32  =  —4ip2q 

(33  ~  —  ip ( r  —  2 qq  )  =  (43  ==  —tu  =  — f«« 
tjs  =  —  2ir<7' 
ti2  =  — 2iTg 
tu  =  -4  ip2q' 
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h\  —  ~tq/02  —  tn  —  ts4 - tit 

£s i  =  -m(T2  +  4 p2qq')  =  £« 

U  2  =  ~4fiTpq 

U 3  =  -m(  r2  -  4p2qq')  =  -£* 

f«  =  ~4mTp<7  ' 

t-n  =  ttt  =  ta  =  hi  =  0 

where 


r  =  2p2  -  1/02. 

After  removal  of  a  common  factor  of  qq',  the  inverse  T„~l  is  seen  to  be  just  a 
rearrangement  of  the  elements  of  Tn 


£«1 

^31 

^31 

^21 

ftl 

0 

0 

‘"^15 

“*33 

“*^33 

^21 

“fl3 

k 3 

—  ^51 

^33 

*33 

^21 

^13 

^€2 

0 

~^42 

““^32 

0 

“^12 

^61 

^31 

^31 

tn 

